A simple reaction-diffusion population model on scale-free networks 
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We study a simple reaction-diffusion population model (proposed by A. Windus and H. J. Jensen, J. Phys. 
A: Math. Theor. 40, 2287 (2007)) on scale-free networks. In the case of fully random diffusion, the network 
topology does not affect the critical death rate, whereas the heterogenous connectivity makes the steady popu- 
lation density and the critical population density small. In the case of modified diffusion, the critical death rate 
and the steady population density are higher, at the meanwhile, the critical population density is lower, which is 
good for survival of species. The results are obtained with a mean-field framework and confirmed by computer 
simulations. 



Recently, Windus and Jensen 01 0] introduced a model for 
population on lattices with diffusion and birth/death accord- 
ing to 2A — > 3 A and A — > <j> for an individual A . They found 
that the model displays a phase transition from an active to an 
absorbing state which is continuous in 1 + 1 dimensions and 
of first-order in higher dimensions [1]. They also investigated 
the importance of fluctuations and that of the population den- 
sity, particularly with respect to Allee effects in regular lattices 
J2]. It was found that there exists a critical population density 
below which the probability of extinction is greatly increased, 
and the probability of survival for small populations can be 
increased by a reduction in the size of the habitat yfl . 

In the study of complex networks [3], an important issue is 
to investigate the effect of their complex topological features 
on dynamical processes |4|], such as the spread of infectious 
diseases [5] and the reaction-diffusion (RD) process @]. For 
most real networks, the connectivity distribution has power- 
law tails P(k) ~ k~" ( , namely, a characteristic value for the 
degrees is absent, hence the scale-free (SF) property. In this 
Brief Report, we shall study the simple RD population model 
|U] on SF networks. 

In an arbitrary finite network which consists of nodes i = 
1, . . . , N and links connecting them. Each node is either oc- 
cupied by a single individual (1) or empty (0). We randomly 
choose a node. If it is occupied, the individual dies with prob- 
ability pd, leaving the node empty. If the individual does not 
die, a nearest neighbor-node is randomly chosen. If the neigh- 
boring node is empty, the particle moves there; otherwise, the 
individual reproduces with probability pb producing a new in- 
dividual on another randomly selected neighboring node, con- 
ditional on that node being empty. A time step is defined as 
the number of network nodes N. In homogeneous networks 
(such as regular lattices and Erdos-Renyi (ER) random net- 
works [7]), the mean-field (MF) equation for the density of 
active nodes p(t) is given by J2l 



dp(t) 
dt 



-p d p(t)+p h {l- Pd )p 2 (t)(l - p{t)), (1) 
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which has three stationary states 



po = 0, p± 



l±Wl- 



4p d 



Fb(l ~Pd) 



(2) 



For 4pd > Pb(l — Pd), Po is me oni y rea l stationary state, 
and one can obtain a critical death rate pd c — Pb/(4 + Pb) 
which separates the active phase representing survival and the 
absorbing state of extinction. Simple analysis shows that p + 
and po are stable stationary states, whereas p_ is unstable and 
therefore represents a critical density p c below which extinc- 
tion will occur in all cases. Thus, forpd < Pd , one can write 

Qi 



pit) 



o 



P(t) < Pc 
P(t) > Pc 



as t 



(3) 



At pd = pd c , the stationary density jumps from 1/2 to 0, re- 
sulting in a first-order phase transition. 

In order to study analytically this process on SF networks 
in which the degree distribution has the form P(k) ~ fc~ 7 
and nodes show large degree fluctuations, we are forced to 
consider the partial densities Pk(t), representing the density 
of individuals in nodes of degree k at time t @]. To ob- 
tain a rate equation for pk(t) , we use a microscopical ap- 
proach which has been applied in diffusion-annihilation JUl 
and multicomponent RD processes on SF networks J^]. Let 
rij(t) be a dichotomous random variable taking values or 
1 whenever node i is empty or occupied by an individual A, 
respectively. Using this formulation, the state of the system 
at time t is completely defined by the state vector n(i) = 
{ni(t), n 2 (i), • • • , njv(t)}. The evolution of n(t) after a time 
step can be expressed as 



rn(t+l) = m{t)r]+ [1 -rii(t)]£, 



(4) 



where r\ and £ are dichotomous random variables, taking val- 
ues of or 1 with certain probabilities p and 1 — p, respectively, 



0; p = p d + (1 - Pd) 



i; l-P, 



1 - i^ay%-(t) 



' (5) 
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0; 1-p 



(1 - Pd)aijrij(t) 
A, ; 



3 / 



(6) 



it occupied. One is that its survival neighbors will move to 
i with probability J2j ^^fc^"^^ anc ' tne other is that its 
survival neighbor j selecting an occupied neighbor I (with the 
probability j- J^i a jini{t))) reproduces a new individual on 



(i-PdKj"j(t) p b 



Obviously, p is the probability that an occupied node i be- *> men we nave me term 
comes empty. If node i is empty, there are two factors causing Taking the average of Eq. ©, we obtain 



(m{t+l)\n(t)) = m(t)(l - p d )l- J2 a ij n j (t) + (l-n,-(t)) 



E (l -Pd)aijrij(t) 
fc, 



(7) 



r 



which describes the average evolution of the system, condi- 
tioned to the knowledge of its state at the previous time step. 
In the MF approximation, (ni(i)rij(t)) = (ni(t))(rij(t)) and obtain 
(n,i(t)n,j(t)ni(t)) = (ni(t))(nj(t)}(ni(t)). Thus, after mul- 



tiplying Eq. <[Vj by the probability to find the system at state 
n at time t and summing for all possible configurations, we 



Pi (t + l) = (n i (t + l)>=p i (i)^ T ^E a ^(*) + ( 1 -^(*)) 



\ - (1 -Pd)aijPj(t) „ , Pb ^ f „ 

i J J 1 



We assume that nodes with the same degree are statistically 
equivalent, i.e., 



Pi (t) ee p k (t) V^ G V(k), 



(9) 



and have 



E a « = E E ^ = j2kp(k'\k) wevw, 

j k' jev(fc') fc' 

(10) 



where V(fc) is the set of nodes of degree fc. 

We split the sum with index j into two sums over fc' and 
V(fc'), respectively. The double sum over ay is related to the 
conditional probability P(k'\k) that a node of given degree 
fc has a neighbor which has degree fc'. In present work, we 
restrict ourselves to the case of uncorrected networks in the 
following, in which the conditional probability takes the sim- 
ple form P(k'\k) = fc'P(fc')/(fc). Thus, from Eq. © and 
after some formal manipulations, we obtain 



p k (t + 1) = p fc (t)(i - Pd)e( P (t)) + jjt(i - Pd)(i - p*(t))^(*)[i +Pb®(p(t))}, 



(ii) 



r 



where p(t) is the total density of active individuals and From Eq. ( fTTT i. we can obtain the following rate equation 
0(p(t)) is the probability that any given link points to an oc- 
cupied node 



Q(p(t)) = ^J2 kp Wp^- 



(12) 



- - Pk (t) + Pk (t)(i- Pd )Q(t) + A(i _ Pd )(i _ p k (t)) P (t)[i+ Ph e(p(t))}. 



(13) 
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FIG. 1: Simulation results of j- — 1 against the reciprocal of k in 
log-log scale on uncorrelated random SF networks with exponent 
7 = 3.0, k min — 3 and TV = 10 3 . The birth rate is p h = 0.5. All 
the plots recover the form predicted in Eq. t!4t . 



Imposing stationarity dtPk(t) = 0, we obtain 



Pk 



-^{l-p A )p{l+p h ®{p)) 



cosh + ^l 



Pd )p{i+ Ph e(p)) 



(14) 



This set of equations imply that the higher the node connectiv- 
ity, the higher the probability to be in an occupied state. This 
inhomogeneity must be taken into account in the computation 
of O(p). Multiplied Eq. ( fT3l by P(k) and summing over k, 
we obtain a rate equation for p(t) 

dp(t) 



dt 



-p(t) Pd + pb(i - Pd)[i - e(p(t))}p(t)e(p(t)), 

(15) 

Notable, the above equation is consistent with Eq. (fl} by im- 
posing that Q(p(t)) = p(t) in homogeneous networks. It also 
has three stationary states 



P S o F - 0, 



e(pT) =p± = 



i 



1±4 1- 



4p d 



The critical death rate pd c 



(16) 

Pb/(4 + Pb) is the same as that 



in homogeneous networks. Imposing, naturally, > 0, 



we find that p+ and p$ are stable stationary states, whereas 
p_ is unstable and therefore represents a critical density F 
below which extinction will occur in all cases. In SF net- 
works, the higher the node connectivity, the higher the prob- 
ability to be in an occupied state (Eq. (TBli). The presence of 
nodes with very large degree results in that p~± < Q(p± F ). 
We conclude that both the population steady state and the crit- 
ical population in SF networks are smaller than those in ho- 
mogeneous networks. 

The numerical simulations performed on uncorrelated SF 
networks confirm the picture extracted from the analytic treat- 
ment. We construct the uncorrelated SF network by the algo- 
rithm [10] with minimum degree fc m ; n = 3 and size N = 



FIG. 2: (a) Stationary population densities p against the death rate 
Pd- (b) For different death rate pa, there are different critical initial 
population densities p(0). All the data were obtained with pb = 0.5 
and TV = 1000. MF results for the homogenous network (line) and 
simulation results for the ER random network (triangles), uncorre- 
lated random SF networks with exponent 7 = 3.0 (circles) and 2.5 
(squares). 



1000. The simulations are carried out on an initially fully oc- 
cupied network for obtaining pk and steady states. To find the 
critical population density, due to its instability, we instead use 
the initial population density p(Q) and find the value of pa that 
separates the active and absorbing states. The prevalence pk is 
computed and averaged over 100 times for each network con- 
figuration, which are performed on 10 different realizations of 
the network. Figure Q] shows the behavior of the probability 
Pk that a node with degree k is occupied. The numerical value 
of the slope in log-log scale is about 0.98, which is in good 
agreement with the theoretical value 1 in Eq. ( [Pil l. In Fig. 
[2] for ER networks with average degree (k) — 14 and size 
N = 1000, both plots of the stationary population density 
p and the critical initial population density p(0) are consis- 
tent with the MF results in homogenous networks (Eq. (0). 
For SF networks, the population and the critical initial popu- 
lation density in steady states are smaller than that of homo- 
geneous networks, which agrees with Eq. ( fTol l. Furthermore, 
the more heterogeneous the SF network, i.e., the smaller de- 
gree exponent of the SF network, the smaller the densities. 
Noting that in despite of the different network topologies, the 
critical death rate pd c is changeless, which is the prediction of 
Eq. <H§). 

In the above model, the neighboring node is chosen with 
full randomness, and we call this fully random diffusion. In 
the following, we shall redesign the diffusion strategies of the 
population model on SF networks. If the randomly chosen 
particle does not die, a nearest neighbor node j is randomly 
chosen with a probability proportional to k". If the neigh- 
boring node j is empty, the particle moves there; otherwise, 
the particle reproduces with probability p\, producing a new 
particle on another neighboring node I which is chosen by a 
probability proportional to fcf , conditional on that node being 
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FIG. 3: (a) Stationary population densities p against the death rate 
Pd- (b) Critical initial population densities p(0) versus the death rate 
p d . All the data were obtained with p h = 0.5 and N = 1000. MF 
results for the homogenous network (line), simulation results in the 
uncorrelated random SF networks with exponent 7 = 3.0 with fully 
random diffusion (triangles), and the modified diffusion under a = 1 
and ft — —2 (circles). 



empty. We can also write the MF equation in the uncorrelated 
random SF network for this modified diffusion case 



dp(t) 
dt 

where 



= -p(t)pd +p b (i -Pd)(i - e p (p{t)))p(t)e a ( P (t)), 

(17) 

Q a (p(t)) = Yl k 1+a P(k) Pk (t). (18) 



In the extinction state p(t) = 0, it is natural that Q a (p(t)) = 
0. We can obtain at lowest order in p(t), <d a (p(t)) ~ Ap(t). 
Similarly Qp(p(t)) ~ Bp(t) can also be obtained at lowest 
order in p(t), where A and B are coefficients. In the steady 
state, following previous analysis, we can easily get the criti- 
cal death rate 



Pdc = 



Ph 



(4|+Pb) : 



(19) 



which can be changed by the ratio B/A, If B < A, the 
modified model has larger critical death rates than the orig- 
inal. Similar to Eq. (fT4l) . the partial density takes the form 



1) 



1/D(k) with D(k) 



xPb®c 



Considering that < pt, and 9 Q < 1, we can negate the sec- 
ond term in D(k). If a > — 1, the higher degree node has 
larger partial density pk. As ft < a, we get B < A; On the 
other hand, if a < — 1, the higher degree node has smaller p}.. 
As ft > a, we have B < A. 



In Fig. [3] the diffusion coefficients are a = 1 and ft = —2. 
From the previous discussion, we obtain B < A. Thus, the 
modified model has a larger critical death rate. Furthermore, 
it has the larger steady population density and the smaller crit- 
ical population density. From the view of conservation, it is 
better that a population system has a larger critical death rate, 
a larger steady population density and a lower critical popu- 
lation density at the same time. Our modified model has this 
nice property under the condition of B < A. 



In summary, we have studied a simple RD population 
model on SF networks by analytical methods and computer 
simulations. We find that in the case of fully random diffusion, 
the network topology can not change the value of the critical 
death rate pd c ■ However, the more heterogenous the network, 
the smaller steady population density and the critical popula- 
tion density. For the modified diffusion strategy, we can ob- 
tain the larger critical death rate and the higher steady popu- 
lation density, at the meanwhile, the lower critical population 
density, which is good for the specie's survival. In present 
work, we consider the population model which has only one 
specie, it may be more interesting to investigate population 
models having multi-species with predator-prey, mutualistic, 
or competitive interactions in complex networks, which is left 
to future work. 
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